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ABSTRACT 

Thin, Keplerian accretion disks generically become gravitationally unstable at large 
radius. I investigate the nonlinear outcome of such instability in cool disks using razor- 
thin, local, numerical models. Cooling, characterized by a constant cooling time r c , 
drives the instability. 

I show analytically that, if the disk can reach a steady state in which heating by 
dissipation of turbulence balances cooling, then the dimensionless angular momentum 
flux density a = ((9/4)7(7 — 1)Qt c )~ . Numerical experiments show that: (1) if r c > 
3S7 1 then the disk reaches a steady, gravito-turbulent state in which Q ~ 1 and cooling 
is balanced by heating due to dissipation of turbulence; (2) if r c < 3f2 -1 , then the 
disk fragments, possibly forming planets or stars; (3) in a steady, gravito-turbulent 
state, surface density structures have a characteristic physical scale ~ 64GS/0 2 that is 
independent of the size of the computational domain. 

Subject headings: accretion, accretion disks, solar system: formation, galaxies: nuclei 



1. Introduction 

It has long been realized that the outer reaches of accretion disks around active galactic nuclei 
(AGN) and young stellar objects (YSO) may be gravitationally unstable (for a review see, for AGN: 
Shlosman et al. (1990); YSOs: Adams &: Lin (1993)). Instability in a Keplerian disk sets in where 
the sound speed c s , the rotation frequency S7, and the surface density £ satisfy 

Q = -"T^yT < Qcrit ~ 1 (1) 

(Toomre 1964; Goldreich & Lynden-Bell 1965). Here Q cr it = 1 for a "razor-thin" (two dimensional) 
fluid disk, and Q cr it = 0.676 for a finite thickness isothermal disk (Goldreich & Lynden-Bell 1965). 
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The instability condition (1) can be rewritten, for a disk with scale height H ~ c s /f2, around a 
central object of mass M*, 

M disfe > -M*, (2) 

where M disk = 7rr 2 £. 

In a steady-state disk whose heating is dominated by interior turbulent dissipation and whose 
evolution is controlled by internal transport of angular momentum, the accretion rate M = 3irac^T,/fl, 
where I have used the a formalism of Shakura & Sunyaev (1973). Then 

M>^ = 7.1xlO- a ( n ^ I ) 3 M sy r- (3) 

implies gravitational instability (e.g. Shlosman et al. (1990)). This does not apply to disks domi- 
nated by external torques (e.g., a magnetohydro dynamic [MHD] wind) or disks heated mainly by 
external illumination. 

For a young, solar mass star accreting from a disk with a = 10 -2 at 10 _6 MQyr _1 equation 
(3) implies that instability occurs where the temperature drops below 17 K. Disks may not be this 
cold if the star is located in a warm molecular cloud where the ambient temperature is greater 
than 17 K, or if the disk is bathed in scattered infrared light from the central star (although there 
is some evidence for such low temperatures in the solar nebula, e.g. Owen et al. (1999)). If the 
effective value of a is small and heating is confined to surface layers, however, as in the layered 
accretion model of Gammie (1996a), then instability can occur at much higher temperatures. 

AGN disk heating is typically dominated by illumination from a central source. The tempera- 
ture then depends on the shape of the disk. If the disk is flat or shadowed, however, and transport 
is dominated by internal torques, one can apply equation (3). For example, in the nucleus of NGC 
4258 (Miyoshi et al. 1995) the accretion rate may be as large as 10 -2 M yr _1 (Lasota et al. 1996; 
Gammie et al. 1999). Equation (3) then implies that instability sets in where T < 10 4 (a/10~ 2 ) K. 
If the disk is illumination-dominated, however, then Q fluctuates with the luminosity of the central 
source. 

The fate of a gravitationally unstable YSO or AGN disk depends on how it arrived in an unsta- 
ble state. To understand why, consider an analogy with convective instability in stellar structure 
theory. The evolution of stellar models with highly unstable radial entropy profiles are of little 
interest, because convection prevents such models from being realized in an astrophysically plausi- 
ble setting. Similarly, highly unstable disks may be irrelevant because the action of the instability 
would prevent one from ever arriving in such a state. Gravitational instability must be "turned 
on" in a natural way. 

An initially stable Keplerian disk can be driven unstable by an increase in surface density (e.g. 
Sellwood & Carlberg (1984)) or by cooling. In an a disk model, the surface density changes on the 
accretion timescale ~ (r/H) 2 (aQ)~ 1 (more rapid variation of surface density could be obtained by 
dumping material onto the disk or by application of a direct magnetic torque). The temperature, 
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by contrast, changes on the cooling time, ~ (a>Q) 1 . This suggests that in cool disks {r/H » 1), 
cooling is the dominant driver of gravitational instability. 

Once gravitational instability sets in, the disk can attempt to regain stability by rearranging 
its mass to reduce E, or by heating itself through dissipation of turbulence. Dissipation can occur 
directly through shocks or indirectly via a turbulent cascade to the viscous scale. If one can model 
the effects of "gravito-turbulence" on a cool disk as an a viscosity, then mass shifting can be accom- 
plished only on the accretion timescale (r/-£f) 2 (af2) _1 . Heating occurs on the thermal timescale 
{aVL)~ l . This suggests that cool disks will attempt to reestablish stability through dissipation of 
turbulence. 

Let us assume that cooling drives the disk toward instability, and that the disk tries to recover 
through turbulent dissipation. There is now the possibility of a feedback loop. If Q is too large 
then heating by gravito-turbulence is weak and the disk cools toward instability. If Q is too small, 
heating by gravito-turbulence is strong and the disk is pushed back toward marginal stability. The 
feedback loop acts as a thermostat to maintain Q ~ 1, just as convection in stars drives the entropy 
profile toward marginal stability. 

The effectiveness of the feedback loop in maintaining Q ~ 1 is controversial. Lin & Pringlc 
(1987, 1990) have considered disk models with Q C 1. Indeed it is not clear that the feedback is 
strong enough to preserve Q ~ 1. Cooling might be so rapid that the disk fragments before it can 
heat (Shlosman k, Begelman 1989), leading to formation of gas giant planets or stars. 

The feedback loop is an old idea in the theory of galactic disks (Goldreich & Lynden-Bell (1965); 
there heating due to star formation balances cooling) and in accretion disk theory (Paczyhski 1978). 
Numerical experiments on cooled collisionless disk models (Anthony & Carlberg 1988; Carlberg & 
Freedman 1985; Tomlcy et al. 1991, 1994) do yield disks with Q ~ 1, and clumping with sufficiently 
strong cooling. Observations of star forming regions of galactic disks also show Q ~ 1 in the gas 
(Kennicutt 1989). 

Numerical experiments with fluid disks often assume a fixed temperature profile or a fixed 
entropy profile (Boss 1997, 1998; Nelson et al. 1998). A fixed temperature profile is relevant to an 
illumination-dominated disk. Other work considers the adiabatic evolution of an initially unstable 
state (Pickett et al. 2000; Laughlin & Rozyczka 1996; Yang et al. 1991). Recent work by Nelson et 
al. (2000) is most comparable to that presented here. It examines the evolution of cooling, global 
models of the disks, but with a more complicated treatment of the cooling function. 

In this paper I consider the nonlinear evolution of cooling, self-gravitating fluid disks via 
numerical experiment. The disk model is razor-thin (two dimensional) and local. In §2 I describe 
the model, then show that its angular momentum flux density can be calculated analytically (§3), 
provided that the disk does not fragment. The results of numerical experiments are described in 
§4. Conclusions are given in §5. 
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The direct, numerical evolution of a global self-gravitating disk model over many dynamical 
times is prohibitively expensive. To develop a converged, calculable, numerical model, and to 
reduce the number of free parameters associated with a global disk model, I make two simplifying 
assumptions: that the disk is cool (i.e. thin, so H/r ~ c s /(Qr) <C 1) and that it razor-thin. 

A cool disk can be modeled using a rigorous expansion of the equations of motion through 
lowest order in c s /(Clr). This model is called the "local model" or "shearing sheet." The expansion 
proceeds as follows. Choose a fiducial point (r ,(p + £l(r )t) that corotates with the disk. Define 
a set of local Cartesian coordinates x,y = r — r a , r (<p — (j) — Q(r )t), and expand the equations of 
motion to first order in |x|/r D . To do this, one needs to assume that the departures from circular 
orbits in the disk |<5v| ~ c s <C Qro, and similarly that the perturbed potential 5(j) ~ c 2 s <C (Qro) 2 . 
With this assumption, the local model describes the evolution of the disk near the fiducial point. 
The self-consistency of these orderings can be tested in the nonlinear evolution. 

I will also use a "razor-thin" disk approximation that treats the disk as two dimensional. There 
is good agreement between razor-thin models and full three dimensional evolution in the linear 
theory of polytropic slender tori (Goldreich et al. 1986), in nonlinear simulations of polytropic tori 
(Hawley 1990), and in nonlinear simulations of gas flow in spiral arms with an isothermal equation 
of state (Tubbs 1980). In general there can be no rigorous mapping between razor-thin disks and 
full three-dimensional systems because the vertical structure contains degrees of freedom that are 
absent in the razor-thin model. 

The equations of motion in the local model read, where v is the velocity, P is the (two- 
dimensional) pressure, and <\> is the gravitational potential with the time-steady axisymmetric 
component removed: 

^ = -X^-2ftx v + 3n 2 xe x - V0. (4) 

For constant pressure and surface density, v = — ^Qxe y is an equilibrium solution to the equations 
of motion. This linear shear flow is the manifestation of differential rotation in the local model. 

The equation of state is 

p = (7 - m (5) 

where P is the two-dimensional pressure and U the two-dimensional internal energy. I will adopt 
7 = 2 throughout. The two dimensional (2D) adiabatic index 7 can be mapped to a 3D adiabatic 
index T in the low-frequency (static) limit. For a nonself-gravitating disk 7 = (3r — 1) /(T + 1) (e.g. 
Goldreich et al. 1986, Ostriker et al. 1992). For a strongly self-gravitating disk, one can show that 
7 = 3- 2/r. 

If the evolution is strictly adiabatic and the initial conditions are isentropic, then p 
with K = const. It follows that the "potential vorticity" 

V x v + 2Q 
* = S 



= KYP 1 
(6) 
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obeys D^/Dt = 0. There is a close connection between the conservation of potential vorticity 
and the stabilizing effects of rotation (Lynden-Bell 1966; Hunter 1964). Processes that cause the 
potential vorticity to evolve can compromise rotational support of the disk at long wavelengths 
(Gammie 1996b). 

The internal energy equation is 

^ + V-(C/v) = -PV-v-- (7) 

at t c 

The final term is the cooling function, and r c is the cooling time. In general r c = r c (S, U, ft). I will 
make the simplifying assumption that r c = const. Notice that there is no heating term; heating is 
due solely to shocks. Numerically, the fluid is heated by artificial viscosity in shocks. 

The gravitational potential is determined by the razor-thin disk Poisson equation: 

V 2 <p = 4itGZ5(z). (8) 
For a single Fourier component of the surface density this has the solution 

= _^*E k e ik - x -IH (9) 

A finite thickness disk has weaker self-gravity, but this does not qualitatively change the dynamics 
of the disk in linear theory (Goldreich & Lynden-Bell 1965). 



2.1. Boundary Conditions, Numerical Methods, and Tests 

I use the "shearing box" boundary conditions, described in detail by Hawley et al. (1995). They 
apply to a rectangular domain of size L x by L y , although I will assume, unless stated otherwise, that 
L x = L y = L. The boundary conditions may be written, for dependent variable / = (X, U, v x , 5v y = 

v y + |rix), 

f(x,y) = f(x,y + L) (10) 
f(x,y) = f(x + L,y-^QLt) (11) 

These boundary conditions have been used to study accretion disks (Hawley et al. 1995), galactic 
disks (Toomre & Kalnajs 1991), and planetary rings (Wisdom &; Tremaine 1988). 

I integrate the governing equations using a self-gravitating hydrodynamics code based on ZEUS 
(Stone & Norman 1992). ZEUS is a time-explicit, operator-split, finite-difference method on a 
staggered mesh. It uses an artificial viscosity to capture shocks. My implementation has been 
tested on standard linear and nonlinear problems, such as sound waves and shock tubes, with 
acceptable results. 
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The transport scheme differs significantly from the basic ZEUS algorithm. I divide the trans- 
port step into two pieces: one due to the mean velocity — ^Qxe y and the other due to the remaining, 
perturbed velocity. The mean velocity transport is done by linear interpolation, while the perturbed 
velocity transport is done as in ZEUS. 1 The "ghost zones" outside the boundaries are likewise set 
using linear interpolation. This algorithm has the advantage that the timestep is limited not by 
the large shear velocity at the radial edges of the grid, — |fi(L/2), but by the perturbed velocities 
only. Longer timesteps can therefore be used, and numerical diffusion is reduced and more nearly 
uniform across the grid. 

I solve the Poisson equation using the Fourier transform method, modified for the shearing 
box boundary conditions. The density is periodic in shearing coordinates x,y', where 

y' = y- hlx{t - t p ) (12) 

and t p = (§fi) _1 NINT(i§fi). Thus when t = t p , which happens at intervals At = (f^)~\ the 
model is periodic in x,y. I linear interpolate £ onto a grid in x,y', solve for the potential via 
Fourier transform, linear interpolate the potential back onto a grid in x, y, and then difference to 
obtain the acceleration. 

There is one subtlety involved in the solution of the Poisson equation. In the Fourier domain, 
the grid of available wavenumbers forms a parallelogram in the k x , k y plane whose shape changes 
with time. I use only those wavenumbers with |k| < (2)~ 1 / 2 (N/2)(2ir/L) in the solution, where 
N x N is the numerical resolution of the model. This restricts the gravitational kernel to the 
largest circular region in Fourier space that is always available. This restriction ensures that the 
gravitational force is approximately isotropic on small scales. 

I have tested my algorithm on three problems. First, I tested the transport piece of the code. I 
evolved, for one rotation period, initial conditions in which v x (t = 0) = const., P = 0, and a square 
of enhanced density was placed in the middle of the grid. The square sheared out, crossed a radial 
boundary, and returned to its initial position. The density enhancement diffused somewhat, but no 
features were introduced at the boundary crossing. Second, I tested the code's ability to conserve 
potential vorticity. I evolved initial conditions containing a sinusoidal velocity perturbation for 
several shear times. Only modest diffusion of potential vorticity was observed. Finally, I checked 
the evolution of a sinusoidal density perturbation against linear theory. 

A comparison of code output with linear theory is shown in Figure 1. The test model has 
L = 40G£o/^ 2 j Q = 1) an d no cooling. The initial conditions are v = — ^ftxe y and S = So + 
<5£cos(k • x), with k x = —2(2tt/L), k y = (2ir/L), and <5X/£o = 5 x 10 -4 . Linear theory (see 
Gammie (1996b), and references therein) gives the evolution of the amplitude 5T,(t) of a shearing 
wave cos(k • x) where k x (t) = k x (0) + ^Qk y t and k y = const.. The heavy solid line in the figure 
shows the linear theory result. The other lines show the evolution of the test models, which start 



X I thank J. Goodman for suggesting this. See also Masset (1999). 
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at N = 32 (poorest agreement with linear theory) and run by factors of 2 up to N = 256. Since the 
radial wavenumber of the disturbance increases with time, the wave eventually becomes wrapped up 
and unresolved. At the end of the run (t = 10ft" 1 ) the radial wavelength is L/13, or 2.46 x (AT/32) 
zones. Evidently the numerical integration converges on the linear theory result. 

3. Analytical Properties of the Model 

Consider the following Jacobidike "integral" for the shearing box: 

r = J d 3 x E5(z) Q« 2 + 1 + ^ + 1^. (13) 

Here (f>T = — |ft 2 x 2 is the tidal expansion of the effective potential about the fiducial point, and 
(j) is the potential of surface density fluctuations. The integral is taken over the entire simulation 
area and from z = — oo to z = oo. Now evaluate dT/dt using the energy equation, the equations 
of motion, the Poisson equation, and the boundary conditions: 

f - \ aL I dS + S§) - 7. J " 'W. <") 

where the surface integral is taken over one of the radial boundaries (the azimuthal boundaries 
do not contribute because they are periodic). To obtain this result, one must realize that terms 
proportional to D(p/Dt vanish when integrated over the boundaries, while terms proportional to 
d<p/dt = D<j>/Dt - (v • V)0 do not. 

Now average equation (14) over t, x, and y, replacing J x dSf by L(J dzf). If the disk can 
settle into a steady state then (dT/dt) = and 

(^W* (15> 

The left hand side is the shear stress, which is proportional to the angular momentum flux density. 

The shear stress can be divided into a hydrodynamic piece and a gravitational piece. The 
hydrodynamic shear stress is obtained via the usual Reynolds decomposition: the fluid velocity is 
is reduced to the sum of a mean steady part — ^Qxe y and a fluctuating part v x e x + 5v y e y . Then 
the associated "Reynolds stress" is 

H xy = T,v x Sv y . (16) 

Notice that the mean radial velocity is of order c s (H/r) 2 , so in the local model, accurate to first 
order in H/r, the mean radial velocity is zero. Thus the local model can be used to calculate a 
shear stress, but a global model is needed to relate the shear stress to an accretion rate. 

The gravitational piece of the shear stress is 



(17) 
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(Lynden-Bell & Kalnajs 1972). The vertical integral arises because the gravitational field outside 
the disk contributes to the shear stress. The vertical integral can be done analytically in the Fourier 
domain, giving the space-averaged gravitational shear stress: 

(Gxy) = £ ffGA ffif k ' 2 - (18) 

The sum is over all Fourier components. 

It is convenient to express the shear stress in terms of an effective a: 

®=^K(Gxy + H xy ). (19) 
Using equation (15), and U = 0^/(7(7 — 1)), I find 

a= (7(7-1)^7-^ . (20) 



For 7 = 2, 

a =lk- <21) 

It is remarkable that it is possible to calculate a analytically for this simple model. Equation 
(20) does not close the system of equations for disk evolution, however, because in general r c is a 
function of the disk temperature. It is the local analog of the result that the surface brightness 
of the disk is independent of the details of the angular momentum transport mechanism. Fixing 
Q would close the system of equations, but Q is not known a priori. This result, however, can 
serve as a check on numerical work. It also explains the empirical relation between cooling rate 
and angular momentum flux reported by Tomley et al. (1994). 



4. Nonlinear Outcome 
4.1. Standard Run 

Consider the evolution of a single "standard" run, with L = 320GS/O 2 and N = 1024. The 
initial velocities are v x = 5v x , v y = — + 5v y , where <5v is a white noise random velocity field of 
subsonic amplitude (the outcome is nearly independent of the details of these perturbations). The 
initial internal energy is such that Q = 1, and the cooling time is r c = 10O -1 . 

The evolution of the kinetic, gravitational, and thermal energy per unit area, normalized to 
, are shown in Figure 2. The evolution of the gravitational and hydrodynamic components 
of q are shown in Figure 3. The a's have been smoothed to make the plot readable. A snapshot of 
the surface density at t = 50O -1 is shown as a color plot in Figure 4; red is high density and blue 
is low density. 
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As the model evolves the small initial velocity fluctuations grow exponentially and develop into 
nonlinear surface density, velocity, and potential fluctuations. Shocks of a predominantly trailing 
orientation develop and heat the gas, while cooling works to reduce the gas entropy. Density 
structures develop which are sheared out by differential rotation and are also predominantly trailing. 
The tendency of velocity and density structures to take on a trailing figure gives rise to a finite a. 

Eventually the thermal energy of the disk settles down to a steady state in which the disk 
is near the point of marginal stability, as hypothesized by Paczyhski (1978). The mean stability 
parameter (Q) = (c s )fi/7rG(£) fluctuates slightly, but averages 2.46 over the last of the run. 

This is larger than the neutrally stable value (Q = 1), but small enough that the disk remains 
sensitive to nonaxisymmetric disturbances. This confirms the hypothesis that the disk maintains 
Q ~ 1. 

The analytic theory of the last section predicts that a = 2/90 ~ 0.022. Averaging over the 
final 1 of the run, I find a = 0.0247. The difference is due to energy nonconservation. Grid 
scale velocity differences are damped by numerical averaging. That fraction of the energy extracted 
from the shear that does not go into shocks and shock heating winds up in a turbulent cascade and 
is lost at the grid scale (one could close the energy equation by introducing a viscosity, but this 
introduces a new set of dynamical complications that will be treated in a separate publication). The 
sense of the difference can be understood by noting that a measures the rate of energy extraction 
from the shear flow. The extracted energy flows by shock heating into thermal energy of the disk, 
and from there to cooling. In a steady state, all must balance. If some of the energy is lost before 
heating the disk, the loss must be compensated for by an increase in the rate of energy extraction. 



4.2. Effect of Cooling Time 

The fates of gravitationally unstable disks can be divided into two classes, depending on the 
cooling time. I find that for short cooling times, r c < 3fi -1 , the model disk fragments. This is easy 
to understand: for very short cooling times pieces of the disk cool and collapse before they have an 
opportunity to collide with one another and reheat the disk (Shlosman & Begelman 1989). 

Fragmentation is illustrated in Figure 5, which shows a snapshot from a run with TV = 512, L = 
80GE/f2 2 ,T c = 2f2 _1 . In this run the disk initially fragments into 6 bound objects, which then 
undergo collisional agglomeration. Ultimately a single bound object forms which is itself a self- 
gravitating accretion disk. This disk has a rotation frequency fi' 3> £1, and so it is able to sustain 
accretion due to repeated gravitational instability because r c > 3fi /_1 . 

When r c > 30 -1 , the outcome is similar to that of the standard run: repeated local cycles of 
cooling, instability, and shock heating. According to the analytic theory of §3, a = ((9/4)f2r c 7(7 — 
1)) _1 in this regime. Figure 6 shows a vs. r c for r c > 3fi -1 . The points are measurements from 
numerical experiments; the solid line is the analytic result. The numerical experiments give slightly 
high values, for reasons discussed in §4.1. 
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4.3. Locality of Angular Momentum Transport 



The preceding discussion is predicated on the idea that gravitational instability can be studied 
in a local model. It has been argued by Balbus & Papaloizou (1999) (hereafter BP) that the 
long-range nature of the gravitational force precludes such a study. Here I argue in return that, 
if certain conditions are fulfilled, it is self-consistent to study transport of angular momentum by 
gravitational instability in the local model. 

The local model is an asymptotic expansion of the governing equations for a thin accretion disk 
to lowest order in e = c s /(roCl) in the neighborhood of a "fiducial point". It assumes that one is 
studying structures on lengthscales ~ ero, that perturbed velocities ~ ero^l, and that the perturbed 
potential 5 4> ~ e 2 r 2 £l 2 . The Poisson equation is solved self-consistently in a WKB expansion (see 
Shu (1970) for higher order corrections). Even if e <C! 1 there are several astrophysically plausible 
ways the local expansion might fail. 

First, structure could exist in the potential on scales large compared to ero. For example, a 
bar might drive a large-scale disturbance in the neighborhood of the fiducial point. Bars cannot 
be generated self-consistently near the fiducial point because, as I will show below, gravitational 
instability only generates structure on scales comparable to ero. The bar must be generated else- 
where in the disk, where e(r) ~ 1. It is therefore a condition for validity of the local model that 
such large-scale external forcing be absent. 

Second, subtler structures may exist in the velocity field, surface density, or pressure on scales 
larger than ero. For example, BP have argued that gravitational instability can extract energy from 
differential rotation, then transform it into waves that propagate over distances large compared to 
ero before dissipating. Such long-range coupling would invalidate a local treatment. While this is 
a live and interesting possibility for all types of disk turbulence, BP use linear theory to argue that 
self-gravitating disks are qualitatively different. 

Global numerical experiments that contain large scale gradients in disk properties could falsify 
or confirm BP's hypothesis regarding wave transport in self-gravitating disks, but they are beyond 
the scope of this paper. Here I will briefly revisit linear theory to show that self-gravitating disks 
are not qualitatively different from nonself-gravitating disks. I will then use the outcome of the 
numerical experiments to show that long-range correlations in surface density, which might be 
expected to develop in the presence of substantial wave transport, are not present. 

Consider a density wave in a razor-thin Keplerian disk. The disk structure varies only on a 
scale r, and e(r) <C 1. The wave oc exp(i( J* dr'k r (r') + im(j) — iujt)), where k r r 3> 1 and m/(k r r) <C 1. 
The WKB dispersion relation is (to — mQ) 2 = v 2 = f2 2 — 27rGS |/c r | + c 2 k 2 . The doppler-shifted 
frequency v 2 has a minimum at \h r \ = k cr = TIGY>q/c 2 . Then the radial energy flux density, referred 
to an inertial frame, is 




(22) 
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(Shu (1970); Goldreich & Tremaine (1979)). This is the full wave energy flux. BP's "anomalous 
flux", by comparison, is the gravitational component of the energy flux measured in a corotating 
frame. Shutting off self-gravity is equivalent to taking k cr — > in eq.(22). Evidently the wave 
energy flux does not change qualitatively in this limit. 

Wave energy fluxes may nonetheless be present. If they are to change disk structure sig- 
nificantly, however, they must be of the same order as the turbulent energy flux FE, W ave = 
(3/2)aScgrO. If I assume that <5£ ~ So, Q ~ 1, and k ~ k cr , and drop factors of order unity, I find 
that \FE,wave/FE,turb\ ~ (m/a) (c s /(rQ)) (for acoustic waves in a nonself-gravitating disk, c s /(rf2) 
is replaced by l/(\k r \r)). Thus for 

(23) 

a r\l 

the wave energy flux is as important as the turbulent energy flux. 

To proceed further one can only consider the plausibility of a large amplitude, high-m wave 
propagating over significant distances in the disk. Here are two arguments against this. First, a 
density wave can only propagate a distance ~ r/m before it turns into an acoustic wave (k r c s > Q). 
In a finite thickness disk this corresponds to a wavelength smaller than a scale height. If the disk 
is stratified three-dimensional effects will modify the wave (e.g. Ogilvie & Lubow (1999)), and 
the wave is likely to steepen, shock, and dissipate. Second, the gravito-turbulent state contains 
fluctuations that emit, scatter, and absorb waves. If scattering and absorption are strong, as they 
are here, coherent signals are destroyed. Under these circumstances it seems unlikely that energy 
will be transmitted over large scales by waves. 

What can the numerical models tell us about the locality of angular momentum transport in 
self-gravitating disks? I have used two methods to assess the locality of structure in the nonlinear 
outcome of my models. In the first analysis, I calculate the dimensionless autocorrelation function 
of the surface density, £ : 

((r) = "I + ^T7^ / d 2 x' E(r + r')S(r'). (24) 

Coherent wave trains would appear as large-scale correlations in the surface density. Figures 7 and 8 
show the autocorrelation function averaged from a series of 5 snapshots at t = (20, 40, 60, 80, 100)S7 — 1 . 
Figure 7 shows the spatial structure of the correlation function from a run with L = 640GS/J7 2 
and iV = 1024. Evidently density correlations are concentrated in a region that is much smaller 
than the size of the model. Figure 8 shows cuts through the correlation function (along the rays 
marked "short axis" and "long axis" in Figure 7) that confirm this quantitatively. 

Also shown in Figure 8 is the autocorrelation function for a run with L = 320GX/$7 2 and 
N = 512 (the same spatial resolution as the larger model). Differences between the smaller and 
larger model result are small and attributable to sampling noise. The correlation function thus 
appears to depend only weakly on L, at least for L > 320G£/f2 2 and r c = lOf^ 1 . This argues that 
surface density structure is locally determined. 
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In a second analysis, I have calculated which Fourier components of the surface density domi- 
nate the gravitational shear stress. Figure 9 shows the quantity 

da G _ 2 f irGk x k y \i: k \ 2 

Here <f) is an angular coordinate in Fourier space. This is the contribution to the gravitational shear 
stress from Fourier components of the surface density in the annulus between k and k + dk. The 
result is calculated from a model with L = 640G£/r2 2 and N = 1024. Fully 90% of the angular 
momentum transport comes from wavenumbers with k > 5(2ir/L). Thus wavelengths significantly 
smaller than the model size dominate the shear stress. 

Figure 9 can be used to estimate how cool a disk must be for the local model to be applicable. 
If the wavenumber k p k of the maximum in dotc/dk is to satisfy k p kr/(2ir) <C 1, then one needs 

S « °- 12 - (26) 

where I have used (27r//c pfc ) pa and (c s ) = 7.7GS/0. YSO disks have H/r ~ 0.1 and 

so are at best marginally described by a local approximation. Some AGN disk models are much 
thinner and thus may be accurately described with a local model. In disks that violate condition 
(26), such as those studied by Laughlin & Rozyczka (1996), global effects can be important. 



4.4. Effect of Resolution 

Most experiments described in this paper have been run at multiple resolution to determine 
convergence. The resolution at which convergence is achieved depends on what is being measured. 
Consider the convergence of a, as measured numerically. Figure 10 shows a, averaged over the 
final 80ft" 1 of the run, in a model with L = 320GS/O 2 and N = 32, 64, ... , 1024. Evidently a 
has converged. Investigation of models of different physical size shows that convergence requires 
resolution of a fixed physical scale which is pa CE/Q, 2 . 



4.5. Length of Run 

The autocorrelation analysis of §4.3 shows that large-scale structures are not established on 
a dynamical timescale. Suppose, however, that gravito-turbulence acts as a viscosity on long- 
wavelength modes. Then there is the possibility of a secular instability analogous to that of viscous, 
self-gravitating disks discussed by Lynden-Bell & Pringle (1974), and also Gor'kavyi et al. (1989); 
Willerding (1992); Schmit & Tscharnuter (1995); Gammie (1996b). These instabilities grow on the 
viscous timescale, ~ (A/ \2ttH)) 2 (aQ)^ 1 . 

As a preliminary test for secular instability, I have integrated two large (L = 160GS/17 2 ), 
high-resolution (N = 512) models for 10 3 fi _1 . One model had r c = 40O -1 , the other r c = 10S7 1 . 
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The energy evolution of the latter is shown in Figure 11. There is no clear trend in any of the 
energy components over the course of the simulation. The r c = run is similar. If secular 

instability is present, it grows only on longer timescales. 

Returning to the analogy with secular instability of viscous, self-gravitating disks, an inte- 
grations over a time tf can be regarded as a test for secular instability over wavelengths A < 
A c = 2-kH yJatfSl. For the integrations described here A c « 30GS/f2 2 . Thus longer integrations, 
although expensive, hold the possibility of interesting results. 

5. Implications 

What are the implications for the structure of disks around YSOs? Although our results are 
not rigorously applicable, since H/r is typically of order 0.1-0.2 in the outer parts of YSO disks, 
they may provide a guide to the relevant physics. 

In models of unilluminated YSO disks (Ruden Sz Pollack 1991; Gammie 1996a) Q approaches 
unity at a radius r ~ 30 AU, where the temperature is « 20 K. Illumination can change Q(r) 
dramatically (Kenyon & Hartmann (1987); for more recent models see, e.g., D'Alessio et al. (1999); 
Bell (1999)). For sufficiently strong illumination, or sufficiently low mass disk, the entire disk is 
gravitationally stable. As the illumination is reduced some radii become gravitationally unstable. 

At gravitationally unstable radii there are two possible outcomes. If the cooling time of the 
disk is sufficiently long (r c > 3f2 -1 ), then gravitationally driven turbulence will simply enhance the 
shear stress above that available from other sources, such as MHD turbulence. This shear stress 
is generated in a region of order across and all the numerical evidence is consistent with 

its being locally determined and treatable as an a viscosity (but only if H/r <C 0.1). Recurrent 
instability and shock heating maintain (Q) 2.4, somewhat larger than the value required for 
axisymmetric stability. 

If, on the other hand, the cooling time in the disk is short (fir c < 3), the outcome is dramatically 
different. Such short cooling times might be found in a disk heated mainly by external illumination 
that is suddenly switched off (for example, in YSO disks, as an FU Orionis outburst turns off). In 
this case the model disk fragments, as suggested by Shlosman & Begelman (1989) in the context 
of AGN disks. The long-term evolution of these systems is uncertain because it is intrinsically 
global; the fragments undergo collisional agglomeration within bands in radius. Clearly one possible 
outcome is the formation of gas giant planets (e.g. Boss (1997) and references therein), but the 
final mass of the object depends on the global evolution. Global models including more realistic 
cooling functions (such as those developed by Nelson et al. (2000)), and external illumination, are 
clearly desirable. 

AGN disk models also generally become gravitationally unstable at large radius (Shlosman & 
Begelman 1989). They typically have H/r <C 0.1 and so the local approximation is more applicable 
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than in YSO disks. Typically at large radius their cooling time is short (equivalently: the a 
required to prevent instability becomes large). The numerical models presented here suggest that 
fragmentation, and perhaps star formation, is likely under these circumstances. 

I am grateful to Jeremy Goodman, Steve Balbus, Ramesh Narayan, and John Papaloizou for 
their comments. This work was supported by NASA grant NAG 52837, NAG 58385, and a grant 
from the University of Illinois Research Board. 
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Fig. 1. — Evolution of the amplitude, in surface density, of a shearing Fourier component. The 
heavy solid line shows the linear theory result. The light lines show measurements from a numerical 
experiment with N = 32 (poorest agreement with linear theory), N = 64, . . . , 256. 
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Fig. 2. — Evolution of the kinetic, gravitational, and thermal energy per unit area, normalized to 
G 2 S 3 /^ 2 , in the standard run, which has L = 320G£/ft 2 , N = 1024, and r c = 10Q -1 . 
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Fig. 3. — Evolution of the gravitational and hydrodynamic pieces of a in the standard run. The 
curves have been boxcar smoothed over an interval 10S1" 1 to make the plot readable. 
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high density. 




Fig. 5. — Map of surface density in a run with r c = 2CI . Black is low density and red is high 
density. The disk has fragmented and formed two bound objects. These objects eventually collide 
and coalesce. 
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Fig. 6. — Time-averaged a for a series of runs with 4 < t c Q < 40 (points) and the analytic scaling 
of §3 (solid line). 
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Fig. 7. — Dimensionless autocorrelation function of the surface density, averaged over a series of 
snapshots from a run with L = 640GS/J7 2 and N = 1024. The lines show the cuts along which the 
correlation function is sampled for Figure 8. Contours are —0.05, —0.025, 0.025, 0.05, 0.075, . . . 0.3. 
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Fig. 8. — The dimensionless autocorrelation function of the surface density, averaged over a series 
of snapshots from a run with L = and N = 1024. The correlation function is sampled 

along the lines shown in Figure 7. The dashed line is for a run with L = 320GS/fi 2 ; the solid line 
is for a run with L = 640GE/S1 2 . 



-25- 




20 40 60 80 100 

k L/(2 tt) 

Fig. 9. — Distribution of the gravitational component of a over |k|, from a run with L = 640GS/fi 2 
and N = 1024. Fully 90% of the shear stress comes from wavelengths A < L/5. 
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Fig. 11. — Evolution of thermal, kinetic, and gravitational energy in a long (lO 3 ^ 1 ) run, with 
r c = lOfi -1 . No clear trend, indicative of a secular instability, is present. 



